echo on
A = [0     2     4     2
     6     8     9     8
     4     0     0     5
     6     2     4     7];
A0=A;
n =4;
pause
P = eye(n);
Q=P;
pause

% scegliamo il pivot

[elm,indm]=max(abs(A(1:n,1:n)))
[elm2,indc]=max(elm);
indr = indm(indc);
pause
% scambiamo le righe di A
A([1 indr],:)=A([indr,1],:)
% scambiamo le colonne di A
A(:,[1 indc])=A(:, [indc,1])

pause     
% scambiamo le righe di P 
P([1 indm],:)=P([indm,1],:)
% scambiamo le colonne di Q
Q(:,[1 indc])=Q(:, [indc,1])

pause
% inizializzo L
L = zeros(n);
% calcolo gli elementi di L :  L(2:4,1)=A(2:4,1)/A(1,1)
L(2:n,1)=A(2:n,1)/A(1,1)
pause
% calcolo la U  riscrivendola in A
% aggiorno la seconda riga
% U(1,:)=A(1,:)
A(2:n,2:n)=A(2:n,2:n)-L(2:n,1)*A(1,2:n);
A(2:n,1)=0
pause
% verifichiamo
P*A0*Q - (L+eye(n))*A
pause
% ripeto il procedimento sull sottomatrice  A(2:n,2:n)

[elm,indm]=max(abs(A(2:n,2)))
[elm2,indc]=max(elm);
indr = indm(indc);

indr = indr+1;
indc=indc+1;
pause
 % scambiamo le righe di A
A([2 indr],:)=A([indr,2],:)
% scambiamo le colonne di A
A(:,[2 indc])=A(:, [indc,2])

pause     
% scambiamo le righe di P 
P([2 indr],:)=P([indr,2],:)
% scambiamo le colonne di Q
Q(:,[2 indc])=Q(:, [indc,2])

% scambiamo le righe di L
L([2 indr],:)=L([indr,2],:)
pause
% calcolo gli elementi di L 
L(3:n,2)=A(3:n,2)/A(2,2)
pause
% calcolo la U  riscrivendola in A
% 
A(3:n,3:n)=A(3:n,3:n)-L(3:n,2)*A(2,3:n);
A(3:n,2)=0
pause
% verifichiamo
P*A0*Q - (L+eye(n))*A
pause
% ripeto il procedimento sull sottomatrice  A(3:n,3:n)

[elm,indm]=max(abs(A(3:n,3:n)))
[elm2,indc]=max(elm);
indr = indm(indc);
indr=indr+2;
indc=indc+2;
pause
% scambiamo le righe di A
A([3 indr],:)=A([indr,3],:)
% scambiamo le colonne di A
A(:,[3 indc])=A(:, [indc,3])

pause     
% scambiamo le righe di P 
P([3 indr],:)=P([indr,3],:)
% scambiamo le colonne di Q
Q(:,[3 indc])=Q(:, [indc,3])


% scambiamo le righe di L
L([3 indr],:)=L([indr,3],:)
pause
% calcolo gli elementi di L 
L(4:n,3)=A(4:n,3)/A(3,3)
pause
% calcolo la U  riscrivendola in A
%
A(4:n,4:n)=A(4:n,4:n)-L(4:n,3)*A(3,4:n);
A(4:n,3)=0
pause
% verifichiamo
P*A0*Q - (L+eye(n))*A

